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Abstract. We study the gravitational clustering of big bang relic neutrinos onto 
existing cold dark matter (CDM) and baryonic structures within the flat ACDM model, 
using both numerical simulations and a semi-analytical linear technique, with the aim 
of understanding the neutrinos' clustering properties for direct detection purposes. In 
a comparative analysis, we find that the linear technique systematically underestimates 
the amount of clustering for a wide range of CDM halo and neutrino masses. 
This invalidates earlier claims of the technique's applicability. We then compute 
the approximate phase space distribution of relic neutrinos in our neighbourhood 
at Earth, and estimate the large scale neutrino density contrasts within the local 
Greisen-Zatsepin-Kuzmin zone. With these findings, we discuss the implications 
of gravitational neutrino clustering for scattering-based detection methods, ranging 
from flux detection via Cavendish- type torsion balances, to target detection using 
accelerator beams and cosmic rays. For emission spectroscopy via resonant annihilation 
of extremely energetic cosmic neutrinos on the relic neutrino background, we give new 
estimates for the expected enhancement in the event rates in the direction of the Virgo 
cluster. 
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1. Introduction 

The standard big bang theory predicts the existence of 10*^^ neutrinos per flavour in 
the visible universe (e.g., PP). This is an enormous abundance unrivalled by any other 
known form of matter, falling second only to the cosmic microwave background (CMB) 
photon. Yet, unlike the CMB photon which boasts its first (serendipitous) detection 
in the 1960s and which has since been observed and its properties measured to a high 
degree of accuracy in a series of airborne/satellite and ground based experiments, the 
relic neutrino continues to be elusive in the laboratory. The chief reason for this is 
of course the feebleness of the weak interaction. The smallness of the neutrino mass 
also makes momentum-transfer-based detection methods highly impractical. At present, 
the only evidence for the relic neutrino comes from inferences from other cosmological 
measurements, such as big bang nucleosynthesis (BEN) and CMB together with large 
scale structure (LSS) data (e.g., Nevertheless, it is difficult to accept that these 
neutrinos will never be detected in a more direct way. 

In order to design possible direct, scattering-based detection methods, a precise 
knowledge of the phase space distribution of relic neutrinos is indispensable. In this 
connection, it is important to note that an oscillation interpretation of the atmospheric 
and solar neutrino data (e.g., implies that at least two of the neutrino mass 
eigenstates are nonrelativistic today. These neutrinos are subject to gravitational 
clustering on existing cold dark matter (CDM) and baryonic structures, possibly causing 
the local neutrino number density to depart from the standard value of tt-i, = ~ 
56 cm~^, and the momentum distribution to deviate from the relativistic Fermi-Dirac 
function. 

In this paper, we develop a method that will allow us to predict the phase space 
distribution of relic neutrinos in our local neighbourhood at Earth (~ 8 kpc from the 
Galactic Centre), as well as in outer space. The method systematically takes into 
account gravitational clustering of relic neutrinos on scales below ~ 5 Mpc, and can be 
applied to the complete range of experimentally and observationally consistent neutrino 
masses. With these predictions, we determine the precise implications of relic neutrino 
clustering for future direct search experiments. To this end, we note that in earlier 
studies of relic neutrino direct detection, the neutrino number density in our local 
neighbourhood is either assumed to be unrealistically large, or simply left as a free 
parameter (e.g., 01 El E!)- With the emergence of the concordance flat ACDM model 
as the cosmological model of choice, today we are in a position to compute the relic 
neutrino phase space distribution within a well deflned cosmological framework, and 
to contemplate again the prospects for their direct detection in a deflnitive way. Our 
studies here will also be useful for such investigations as relic neutrino absorption [HIHlin] 
and emission jTHl HH 1121 UHl IHI spectroscopy. 

The standard procedure for any gravitational clustering investigation is to solve 
the (1 + 3 + 3)-dimensional Vlasov, or coUisionless Boltzmann, equation using A^-body 
techniques (e.g., [IS Uni El CHI)- However, these techniques are computationally very 
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expensive and necessarily come with limited resolutions. In the context of the cold+hot 
dark matter (CHDM) model, earlier A^-body studies involving neutrinos probe their 
kinematic effects on structure formation from cluster and galaxy abundances on large 
scales (e.g., [HI 1201 1211 123 ), to halo properties on small scales |23I. While the CHDM 
model has fallen out of favour in recent years (see, however, |21]), it is instructive to 
note that the halo simulation of [22] has a formal resolution of only ~ 100 kpc. This is 
clearly inadequate for our considerations, where the scale of interest is of order 1 kpc. 

In the context of the fiat ACDM model, Singh and Ma (hereafter, SM) presented a 
novel approximate method to probe the accretion of neutrinos onto CDM halos at scales 
below ~ 50 kpc |2S1- The salient feature of this study is their use of parametric halo 
density profiles from high resolution, pure ACDM simulations as an external input, while 
the neutrino component is treated as a small perturbation whose clustering depends 
on the CDM halo profile, but is too small to affect it in return. Implementation of 
this approximation requires the neutrino mass density Pi, to be much smaller than its 
CDM counterpart pm- On cosmological scales, we know now from LSS data that the 
ratio Pu/Pm = ^u/^m is at most ~ 0.2 [2]. On cluster/galactic scales, neutrino free- 
streaming ensures that p^/ pm always remains smaller than its cosmological counterpart 
|23j . Thus the approximation scheme, so far, is sound. Furthermore, in order to track 
the neutrino density fluctuations in the most effortless way, SM employed the linearised 
Vlasov equation instead of its full version. Unfortunately, linear methods are known to 
break down when the density fluctuations reach the order of unity. Indeed, in their two 
trial runs with CHDM parameters, the linear results of SM compare favourably with 
iV-body results of [2S1 in the outer part of the halo, where the neutrino overdensity is 
relatively low. The denser inner parts ( < 1 Mpc), however, show marked disagreement. 
This discrepancy renders SM's claim that their complete prescription is able to probe 
neutrino clustering on sub-galactic scales doubtful. 

In the present investigation, we adopt one of the more attractive features of SM's 
study, namely, the use of parametric halo profiles as an external input. However, we 
improve upon their analysis by solving the Vlasov equation in its (almost) full glory 
utilising a restricted, A^-l-body (pronounced: EN-ONE-BODY) method based on the 
following observation: In the limit p^ -C Pm and the CDM contribution dominates the 
total gravitational potential, not only will the CDM halo be gravitationally blind to the 
neutrinos, the neutrinos themselves will also have negligible gravitational interaction 
with each other. This allows us to track them one particle at a time in N independent 
simulations, instead of following particles simultaneously in one single run, as in 
a conventional A^-body study. An obvious advantage of our A^-l-body technique is 
that it requires virtually no computing power when compared with a full scale A^-body 
simulation with the same, large A^ (> 10^). It is also less time-consuming since we 
have done away with the need for a gravity solver (the core of all A^-body techniques). 
In addition, restricted methods such as ours do not suffer from spurious two body 
relaxation, and hence do not require the introduction of an artificial softening length that 
is mandatory in conventional A^-body studies. Lastly, we note that restricted methods 
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have been used extensively in the studies of galaxy interactions (e.g., |2ni), and, when 
properly motivated, should not be seen as inferior to full scale iV-body techniques. 

As a closing remark, let us stress again that our purpose here is not to investigate 
the effects of neutrino mass on cosmology, but rather to address some simple questions 
such as how many relic neutrinos can we realistically expect to find in this very space 
we occupy, what kind of energies do they have, where in the universe can we expect to 
find the highest concentration of relic neutrinos, etc., given what we know today about 
cosmology. In this regard, the analysis we present here is most exhaustive. 

The paper is organised as follows. We begin in section El with an assessment of 
the current observational constraints on the relic neutrino background. In section El 
we introduce the Vlasov equation which is used to track the phase space distribution 
of the neutrinos. Section HI contains a brief discussion of the halo density profiles to be 
employed in our calculations. In section we solve the Vlasov equation for the halo 
models using our improved A^-l-body method for a variety of halo and neutrino masses. 
In section IHl we compute for the same halo models and neutrino masses the neutrino 
overdensities using the linear method of SM and examine its validity. Section [71 deals 
exclusively with relic neutrinos in the Milky Way, and in particular their phase space 
distribution in our immediate vicinity. We discuss in section |S1 the implications of our 
findings for scattering-based detection methods, and we conclude in section El 

2. Observational constraints on the relic neutrino background 

Taking as our basis (i) the fiat ACDM model with (f2m,o,^A,o) ~ (0.3,0.7) and Hubble 
parameter /i ~ 0.7, (ii) neutrino mass splittings inferred from the solar and atmospheric 
data, (Arrig^j^, Am^^-jjj) ~ (10^^,10^'^) eV^, and (iii) the invisible Z decay width from 
LEP which constrains the number of SU{2) doublet neutrinos to three [2Z|, a minimal 
theory of neutrino clustering is fixed only by the absolute masses of the neutrinos m,^. 
The current laboratory limit from tritium j3 decay experiments is mj, < 2.2 eV {2a) 
|28| l^. and should improve to ~ 0.35 eV with the upcoming KATRIN experiment |30j . 
Cosmology also provides a constraint on m^. For three degenerate species, an upper 
bound of Y^rUi, < 1.7 eV {2a) jSH E21 EH] has been inferred from a combined analysis 
of the CMB anisotropy from WMAP jSl] and galaxy clustering from SDSS (SDSS-gal) 

(or from 2dFGRS ISEj), together with an HST prior on the Hubble parameter p?7j . 
(Reference [221 uses also SNIa ^|).| Adding to the fit galaxy bias ^\ and Lya forest 
analyses can tighten the constraint to ^2^^ < 0.42 eV ^ (see also jUl US ESI El!), 
although the robustness of these additional inputs is still contentious. Weak lensing of 
galaxies or of the CMB (1^1 will provide an alternative probe for the cosmological 
implications of massive neutrinos. 

While constraints from cosmology are interesting in their own right, they are also 
highly model dependent, and degeneracies abound. For instance, if tensors and running 

I The mass splittings inferred from the solar and atmospheric neutrino oscillation experiments imply 
that the three mass eigenstates are quasi-degenerate when nii, ^/Afnl^. 
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of the scalar spectral index are allowed, the last J2^u bound relaxes to 0.66 eV |4(Jj . 
Another possibility is an interplay between m^, and N^,, where N^, is the effective number 
of thermalised fermionic degrees of freedom present in the radiation-dominated era, such 
that increasing N^, actually weakens the bound on nii, [33 . For example, a N^, = 6 
model receives a CMB+LSS+priors constraint of (i) Yl,^v < 2.7 eV, if all six particles 
are equally massive, (ii) J2^u < 2.1 eV, if three are massive and the others exactly 
massless, and (iii) J2 < 4.13 eV, if only one is massive |17j. Currently, 1.4 < N^, < 6.8 
is allowed by CMB+LSS+priors [SSI 112 UHl 113 ED! • Future CMB experiments such as 
Planck will be sensitive to AA^^ ~ 0.2 |HniH^ 

Lastly, we note that BBN prefers 1.84 < < 4.54 {2a), in the absence of a 
chemical potential j^S] (see also [301 154j ) . Allowing for a nonzero weakens the 
bounds to 1.3 < A",^ < 7.1 for —0.1 < (^e < 0.3 jHHl- There is no lack of candidates in the 
literature for these extra N^, — 3 degrees of freedom. We shall not list them here. What 
is certain, however, is that they cannot take the form of very large chemical potentials 
in the i/^^t sectors, since large neutrino mixing inferred from oscillation experiments 
ensures that ~ ~ < 0.3, too small to be a significant source of N^, jSHl EZl EH] ■ 

In the present analysis, we assume the neutrinos to constitute exactly three 
thermalised fermionic degrees of freedom, and adopt a conservative mass bound of 

m,<0.6eV, (2.1) 

corresponding to the Ni, = 3 constraint from the WMAP+SDSS galaxy cluster analysis 
|31j . Alternatively, ()2.1|) may be interpreted as a restrictive bound for models with 
extra, non-neutrino relativistic particles (A^^, > 3), or with a significant running spectral 
index, as discussed earlier. 



3. Vlasov equation 

A system consisting of several types of weakly interacting, self-gravitating particles [e.g., 
CDM plus neutrinos] may be modelled as a mult i- component collisionless gas whose 
phase space distributions /j(a;,p, r) obey the Vlasov equation (e.g., [UlEni), 

^ = — + i- — +■• — = (3 1) 

Dt Ot dx dp 

The single-particle phase density fi{x,p,T) is defined so that dNi = fi d^x d^p is the 
number of i particles in an infinitesimal phase space volume element. The variables 

x = r/a{t), p = arUiX, dr = dt/a{t), (3.2) 

are the co moving distance, its associated conjugate momentum, and the conformal time 
respectively, with a as the scale factor and rui the mass of the ith particle species. 
All temporal and spatial derivatives are taken with respect to comoving coordinates, 
i.e., ' = d/dr, V = d/dx.% In the nonrelativistic, Newtonian limit, equation p.lj) is 

§ Unless otherwise indicated, we shall be using comoving spatial and temporal quantities throughout 
the present work. Masses and densities, however, are always physical. 
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equivalent to 

^f^ , P dfi dfi 

— H amiVcj) ■ — = 0, (3.3) 

OT ami ox op 

with the Poisson equation 

VV = 47rGa2^p,,(r)5,(a3,r), (3.4) 

i 

5i{x,T) = P'^f'^^ - 1, pi{x,T) = ^ f d^p fi{x,p,T), (3.5) 
P,{t) a-" J 

relating the peculiar gravitational potential (j){x,T) to the density fluctuations 5i{x,T) 
with respect to the physical mean Pi(r). 

The Vlasov equation expresses conservation of phase space density fi along each 
characteristic {x{t),p{t)} given by 

— = , — = -aniiVcp. (3.6) 

dr ami dr 

The complete set of characteristics coming through every point in phase space is thus 
exactly equivalent to equation 1)3.11) . It is generally not possible to follow the whole set 
of characteristics, but the evolution of the system can still be traced, to some extent, if 
we follow a sufficiently large but still manageable sample selected from the initial phase 
space distribution. This forms the basis of particle-based solution methods. 



4. Halo density profiles and other preliminary concerns 

A "first principles" approach to neutrino clustering requires the simultaneous solution 
of the Vlasov equation ()3.1|) [or, equivalently, the equations for the characteristics ()3.6|) ] 
for both the CDM and the neutrino components. In our treatment, however, we assume 
only the CDM component pm contributes to in the Poisson equation ()3.4|) . and pm to be 
completely specified by halo density profiles from high resolution ACDM simulations. 
We provide in this section further justifications for this approach, as well as a brief 
discussion on the properties of the halo density profiles to be used in our analysis. 

It is well known that after they decouple from the cosmic plasma at T ~ 1 MeV, 
light neutrinos {mi, -C 1 keV) have too much thermal velocity to cluster on small 
scales via gravitational instability in the early stages of structure formation. Accretion 
onto CDM protoclusters becomes possible only after the neutrino velocity has dropped 
below the velocity dispersion of the protoclusters. The mean velocity of the unperturbed 
neutrino distribution has a time dependence of 

(v) ^ 1.6 X 10^ (1 + z) (—] km (4.1) 

where z is the redshift. A typical galaxy cluster has a velocity dispersion of about 
1000 km s~^ today; a typical galaxy, about 200 km s~^. Thus, for sub-eV neutrinos, 
clustering on small scales can only have been a z < 2 event. 
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On the other hand, based on a systematic A^-body study of halo formation in 
a variety of hierarchical clustering cosmologies, Navarro, Frenk and White (hereafter, 
NFW) argued in 1996 that density profiles of CDM halos conform to a universal shape, 
generally independently of the halo mass, the cosmological parameters and the initial 
conditions jnOllS]- This so-called NFW profile has a two-parameter functional form 

Phalo(r) = , I ..^\ (4.2) 

where r<j is a characteristic inner radius, and ps = 4p(rs) a corresponding inner density. 
These parameters and ps are determined by the halo's virial mass Mvir and a 
dimensionless concentration parameter c defined as 

c^"—, (4.3) 

rs 

where rvir is the virial radius, within which lies Mvir of matter with an average density 
equal to Avir times the mean matter density pm at that redshift, i.e., 

Mvir = -^^vhPmO' ?"vir ~ AvirPm.O'^vir 



3^3 



Anpstt r 



ln(l + c) - 



(4.4) 



1 + cJ ' 

where pm,o is the present day mean matter density. The factor Avir is usually taken to 
be the overdensity predicted by the dissipationless spherical top-hat collapse model (5th, 
which takes on a value of ~ 178 for an Einstein-de Sitter cosmology, while the currently 
favoured ACDM model has (5th — 337 aX z = 0.|| 

Within this framework, any halo density profile is completely specified by its virial 
mass and concentration via equations ()4.2|) to ()4.4p . Indeed, the NFW profile in its two- 
parameter form generally gives, for quiet isolated halos, a fit accurate to ~ 10% in the 
range of radii r = 0.01 1 Tvir PH]- Furthermore, NFW argued for a tight correlation 
between Mvir and c, such that the mass distribution within a halo is effectively fixed 
by the halo's virial mass alone. Later studies support, to some extent, this conclusion 
(e.g., inni IMl); halo concentration correlates with its mass, albeit with a significant 
scatter. The analysis of [HI] of ~ 5000 halos in the mass range 10^^ 10^^ M© reveals 
a trend (at 2 = 0) described by 

with a la spread about the median of A(logc) = 0.18 for fixed Mvir. In addition, for a 
fixed virial mass, the median concentration parameter exhibits a redshift dependence of 

c[z) ~ (4.6) 
^ ' l + z 

between 2; = and z = A. 

In their analysis, SM interpreted the set of equations (j4.2p to (j4.6p as a complete 

description of an individual halo's evolution in time. While we do not completely 

II In the original work of NFW [201 E], ''vir is taken to be the radius r2oo, within which the average 
density is 200 times the critical density of the universe, irrespective of the cosmological model in hand. 
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disagree with this interpretation, it should be remembered that the relations (j4.5p and 
fl4.6|) refer only to the behaviour of the statistical mean for fixed virial masses, and do 
not actually describe how individual halos accrete mass over time. A better motivation 
for the equations' use comes from the observation that the physical densities of the inner 
regions of individual isolated halos tend to remain very stable over time (z ~ 2 — > 0) 
pill. This behaviour, as it turns out, can be more or less reproduced by equations ()4.2p 
to ()4.6|) . if we interpret M^^^ as the halo's virial mass today. Merging subhalos tend 
to affect the main halo's density profile only in its outer region. Therefore, provided 
that neutrino clustering becomes important only after z ~ 2 and the halos have had no 
major mergers since then, we are justified to either use these equations here, or simply 
model the CDM halo as a static object in physical coordinates, as a first approximation. 

Finally, we note that there exists in the literature a number of other halo profiles 
(e.g., inni EZI), which, in some cases, provide better fits to simulations than does the 
NFW profile ()4.2|) . However, these profiles generally differ from ()4.2|) by less than 10% 
[IH] so we do not consider them in our study. 



5. A^-l-body simulations 

Using the NFW halo density profile ()4.2p as an input, we find solutions to the Vlasov 
equation in the limit pi, -C Pm by solving the equations for the characteristics ()3.6|) . We 
discuss below the basic set-up. Technical details can be found in the Appendix. 



5.1. Basic set-up and assumptions 

We model the CDM distribution as follows. We assume that throughout space is a 
uniform distribution of CDM. On top of it, sits a spherical NFW halo at the origin. 
In order that the halo overdensity merges smoothly into the background density, we 
extend the NFW profile to beyond the virial radius. Furthermore, for convenience, we 
treat the NFW profile as a perturbation pm{T)6m{x,T), rather than a physical density. 
This simplification should make very little difference to the final results, since the halo 
density is always much larger than the background density. The halo's properties and 
its evolution in time are contained in the set of equations ()4.2|) to ()4.6|1 . For the factor 
Avir, we take a time-independent value Avir = 200, following and SM. We choose 
this somewhat uncommon definition so as to facilitate direct comparisons between the 
results of SM and those of the present study. However, since the choice of Avir affects the 
profile only through [see equation ()A.10|) in the Appendix], we can see immediately 
that using instead the more common Avir = ^th, where 

187r^ + 82y - 39y^ n f ^ ^ ^^ 

St-r^ 7^-Y~\ ' y = ^m{z)-l, (5.1) 

with VLm{z) = Qm,o/{^mfl + ^Afid^) wiU make little difference to the outcome. 

A more important issue is the role played by other CDM structures (i.e., other 
halos and voids) that should realistically be in the surroundings. In our present scheme. 
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these have all been lumped into one uniform background so as to preserve the spherical 
symmetry of the problem. In reality, these structures will induce tidal forces and distort 
the symmetry. However, we expect tidal forces to be important only for clustering in 
the outer part of the halo where the gravitational potential is low and several halos may 
compete for the same neutrinos. Clustering in the deep potential well of the inner part, 
on the other hand, should not be seriously affected. 

For the neutrinos, we take the initial distribution to be the homogeneous and 
isotropic Fermi-Dirac distribution with no chemical potential. 

Up) = T- — Vtt^, (5.2) 

where Tj^^ = 1.676 x 10~^ eV is the present day neutrino temperature. In principle, 
the chemical potential need not be exactly zero. In fact, a positive chemical potential 
(i, should improve, to some extent, the clustering of neutrinos (as opposed to anti- 
neutrinos) by providing more low velocity specimens that cluster more efficiently than do 
their high velocity counterparts. However, this enhancement is necessarily accompanied 
by a suppression of clustering in the anti-neutrino sector, for which a negative Cp = —(u 
tends to deplete the low velocity states. Currently, the upper bound on (i, is 0.3, too 
small to warrant a detailed investigation into a possible "clustering asymmetry" . 

We simulate initial momentum values in the range 0.01 < p/T^^ < 13, which 
accounts for more than 99.9% of the distribution ()5.2|1 . The initial spatial positions of 
the neutrinos range from r = 0, to as far as it takes for the fastest particles to land 
within a distance of 10 h~^Mpc from the halo centre at 2; = 0. We consider a sample 
of three neutrino masses, m^, = 0.15,0.3,0.6 eV, consistent with the bound ()2.H) . and 
a range of halo viral masses, Mvir = lO^^M© — lO^^M©, corresponding to halos of the 
galaxy to the galaxy cluster variety. All simulations model a flat ACDM cosmology, 
with the parameters (fi^.O; ^a,0! ^) = (0.3,0.7,0.7) 

The initial spatial and momentum distribution is divided into small chunks that 
move under the external potential of the CDM halo, but independently of each other. A 
low resolution run is flrst carried out for each set of {m^, Mvir}. AH chunks that end up 
at z = inside a sphere of radius 10 h~^Mpc centred on the halo are traced back to their 
origin, subdivided into smaller chunks, and then re-simulated. The process is repeated 
until the inner ~ 10 h~^kpc is resolved. The initial redshift is taken to be 2; = 3. 
This should be sufficient, since clustering is not expected to be fully under way until 
2; ~ 2. The final neutrino density distribution is constructed from the set of discrete 
particles via a kernel density estimation method [HHl [TOj outlined in the Appendix, with 
a maximum smoothing length of ~ 2 /i^^kpc in the inner ~ 50 /i^^kpc of the halo. 

5.2. Results and discussions 

The basic results of our A^-l-body simulations are presented in Figure ^ which shows 
the neutrino overdensities Uy/fii, for various sets of {ruy, M^i-^}. A companion figure, 
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Figure 1. Relic neutrino number density per flavour, = rii?, normalised to 
TT-:^ — — 56 cm~^, for neutrino masses m,y — 0.6,0.3,0.15 eV and halo virial 
masses indicated in the figure. Results from A^-l-body simulations are denoted by 
red (solid) lines. Dotted lines correspond to overdensities calculated with the linear 
approximation. 
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Figure 2. Mass density ratio pv/ Pm normalised to the background mean pv/pm 
obtained from iV-l-body simulations for various neutrino and halo masses indicated in 
the figure. 



Figure shows the same results expressed in terms of the mass density ratio p^/ pm 
normahsed to the background mean py/pm- 

The essential features of the curves in Figures and |21 can be understood in terms 
of neutrino free-streaming, which causes the n^jny curves to flatten out at small radii, 
and the mass density ratio p^j pm to drop substantially below the background mean. 
(The latter feature also provides a justiflcation for our A^-l-body method.) Both Uyjuv 
and pyj Pm approach their respective cosmic mean of 1 and pyj pm at large radii. Similar 
behaviours have also been observed in the CHDM simulations of reference [22] ■ 
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Figure 3. Dependence of the neutrino overdensity on the halo virial mass for neutrino 
masses nii, — 0.6,0.3,0.15 eV. The circles represent overdensities at 1 h^^kpc, 
pentagons at 10 h~^kpc, squares at 100 h~^kpc and triangles at 1000 h~^kpc. The 
straight lines (oc Mvir) are provided to guide the eye and are not meant to be best fits. 



Naturally, clustering improves with increasing neutrino and/or halo masses. In 
Figure we plot the neutrino overdensities at 1, 10, 100, 1000 h~^kpc for various 
neutrino masses as a function of the halo virial mass M^^^. A similar plot is constructed 
in Figure m with the neutrino mass m^, as the independent variable. In the former, the 
quantity n^/uy — 1 is seen to be roughly proportional to Mvir for a fixed radius and a 
fixed neutrino mass. The Uy/fiy — 1 versus m,^ trend in Figure |3] is more difficult to 
quantify. A roughly dependence can be discerned for some fixed halo masses at 
some fixed radii. Other combinations, however, display noticeably different behaviours. 
Whatever these dependences are, it is interesting to note that they are never shallower 
than oc ruy, or steeper than oc m^. As we shall see later, an ml and an rnl dependence 
for the overdensities can both be motivated from theory. The former can be derived 
from the linearised Vlasov equation (cf. section IHl), while the latter follows naturally 
from phase space considerations (or the so-called Tremaine-Gunn bound, cf. section [7j). 

Lastly, in order to test the robustness of our results, we (i) push the initial redshift 
of the simulations back to 2; = 5, (ii) vary the cosmological parameters within their 
allowed ranges, and (iii) alter the time dependences of some of the halo parameters. 
In all cases, we find, as expected, the heavier masses {m^^, Mvir} to suffer more from 
these variations. For our heaviest set, {my = 0.6 eV, Mvir = 0.7 x 10^^ H'^Mq}, the 
neutrino overdensity changes by about 10 to 20% at small radii, and some 50 to 100% at 




Figure 4. Dependence of the neutrino overdensity on the neutrino mass for various 
halo masses indicated on the plots. The circles represent overdensities at 1 h'^kpc, 
pentagons at 10 h^^kpc, squares at 100 h^^kpc and triangles at 1000 h^^kpc. The 
solid lines correspond to an dependence, dashed lines an dependence, and 
dotted lines an dependence. These lines are provided to guide the eye and are not 
meant to be best fits. 
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r > 5 /i^^Mpc. The gap narrows with smaller neutrino and halo masses. For galaxy size 
halos (Mvir ~ IO^^Mq), the density variations with respect to (i), (ii) and (iii) are no 
more than ~ 10% everywhere. Thus our simulation results are generally quite robust. 



6. Linear approximation 

The linear approximation is often used in the literature to find approximate solutions 
to the Vlasov equation. In the nonrelativistic limit, the pioneering work of Gilbert |Zl] 
has, over the years, been apphed to the study of the pure hot dark matter (HDM) 
model [72], as well as in the analysis of HDM accretion onto nonadiabatic seeds such 
as cosmic strings [ZSl HH IZS| and onto CDM halos j2S] • The procedure consists of first 
switching to a new time variable s = J a~^dT = J a~^dt, and then Fourier transforming 
the a;-dependent functions, 

f{k,p,s) = T[f{x,p,s)], (f){k,s) = J^[(j){x,s)], (6.1) 

to obtain a new differential equation in Fourier space, 

^ + 'JllPf _ ,rn,a\k * Vj) = 0, (6.2) 

(yS TTXy 

where k (/)-k'Vpf = J d^k' k' 0(A;')- Vp/(fc — A;') is the convolution product. The equation 
is said to be linearised when one makes the replacement 

Vp/>)^Vp/o(p)5(fc), (6.3) 

where 6{k) is the Dirac delta function, so that the convolution product becomes a 
simple scalar product kcp ■ Vp/o. Replacement ()6.3|1 is valid as long as the condition 
|Vp(/ — /o)| ^ |Vp/o| holds. In practice, however, the quantity Vp/ is somewhat 
cumbersome to compute, so the "rule of thumb" regarding the linear approximation 
is to abandon it as soon as the spatial density fluctuation S,,{x,s), deflned in (j3.5p . 
exceeds order unity, as emphasised in [SHI d] • Nonetheless, the linear theory has been 
time and again used beyond this putative limit. We shall also apply it to our case, to 
gain physical insight as well as to see how it compares with A^-l-body simulations. 
Equation ()6.2|1 together with the replacement ()6.3j) has a very simple solution, 

/(fc,p,s) = /(fc,p,s,)e-'=--(^-^') 

+ im^k- Vp/o /' ds'a\s)4){k, s')e-'''-''^'-''\ (6.4) 

J Si 

where u = p/niu, Si is some initial time, and we take the initial phase space distribution 
to be isotropic and homogeneous in space, i.e., f{k,p,Si) = 6{k)fo{p). The neutrino 
number density per Fourier mode relative to the mean density is obtained by integrating 
()6.4|) over momenta p, 

h^{k,s) a-^ J d^p f{k,p,s) _ 1 



n^{s) a-^Jd^pfo{p) u^^q 

= 6{k) - e r ds'a\s')^{k, s'){s - s')F 

J Si 



j d^p /(fc,p,s) 

k{s 



(6.5) 
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Table 1. Some distribution functions fo{p) and their corresponding F{q) [equation 
()6.6|l ] that have appeared in the Uterature. The series solution for F{q) for the 
relativistic Fermi-Dirac (FD) function was first derived in A Maxwell-Boltzmann 
(MB) type distribution was adopted in jj^]- The last family of distributions, 
characterised by F{q)^s exponential form, appears in [HHI and |7f)j . 



Distribution 


Mp) 


F{q) 


Relativistic FD 


[l + exp(p/T,,o)]-i 




MB 


exp(-p/r,,,o) 




7 distribution 




exp(-79T^,o) 



F{q)^— I d'pe-'^y,{p). (6.6) 



with 



The correct form for /o(p) should be the relativistic Fermi-Dirac function (|5.2|) . which 
gives for F{q) a series representation [Tlj 

where C(3) — 1.202 is the Riemann zeta function. However, in order to simplify 
calculations and/or to gain physical insight, other forms of fo{p) have also appeared 
in the literature. Some are listed in Tabled along with their corresponding F{q). 
The solution to the Poisson equation (j3.4p in Fourier space is 

m s) = = — . (6.8) 

Substituting this into ()6.5|1 and using the definition 6i^{k, s) = h^{k, s)/ni^{s) — 6{k), we 
obtain for the neutrino density fluctuations 

^ kjs-s'Y 
my 

This is the "master equation" for the linear approach. We solve equation ()6.9j] 
numerically for a variety of neutrino and halo masses. The results are presented in 
Figure HI alongside their A^-l-body counterparts. 

6.1. Further approximations and analytical insights 

Before comparing the two approaches, let us first study the linear approximation for its 
own sake. Consider the master equation ()6.9j) . In the limit F{q) grows much faster than 
a(s) and 5m{k,s), i.e., 

kTyfl Ida 1 dSm ^ 



6y{k, s) ~ 47rG'p„,o / ds'a{s')6^{k, s'){s - s')F 

J Si 



(6.9) 
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equation (|6.9|) may be solved by asymptotic expansion. The resulting approximate 
solution looks somewhat messy at first sight, 

2 



m,. 



kT, 



u,0 . 



3C(3) 



[\n{2)a{s)6Uk,s)- 



n+1 



n 



n=l 



(6.11) 



S - Si 



but may be rendered into a physically transparent form if we first cross out the second 
term, which is well justified since the initial a and 6m should always be much smaller 
than the final ones, and then rewrite the expression as 

kl(s) 



Suik,s) ~ 



-6m(yk, s). 



(6.12) 



Here, kfs is but the free-streaming wave vector, defined as 



k{Js) = 



\ 



4:TTGa{s)p, 



m,0 



r2 



\ 



ATTGa'^{s)prnis) 

cl(s) 



1.5 Ja{s)a 



m,0 



and we identify 



cAs = 



u,0 



3C(3) 



Cu,0 



(6.13) 



m;^a(s) ^ 2 ln(2) a{s) 



km s 



-1 



(6.14) 



a{s) \m^J 

as the neutrino's characteristic thermal speed. 

The functional form of equation ()6.12|1 already tells us something very interesting; 
large Fourier modes in the neutrino density fluctuations are suppressed by a factor 
proportional to relative to their CDM counterparts. This is clearly a manifestation 
of free-streaming, which is responsible for inhibiting the growth of structures on scales 
below Afs = 2n/kis- Furthermore, 6i, has an ml dependence through kfs, meaning that, 
at small scales, a neutrino twice as heavy as another is able to cluster four times more 
efficiently. This ml dependence is reflected, approximately, by both our linear and 
A^-l-body results in Figured and is particularly pronounced at small radii. 

That equation ()6.12j) is a solution of ()6.9j) is contingent upon the satisfaction of 
the condition ()6.10p . which requires, for a fixed neutrino mass, k to be larger than 
some nominal fcmin determined by the rates of change of the scale factor a and of the 
CDM perturbations 6m{k,s). The rate of change of a is a simple and well defined 
function of the cosmological model. The growth rate of Sm,{k,s), on the other hand, is 
usually more complicated. However, because our halos are practically static in physical 
coordinates, this rate (in comoving Fourier space) can only be at most of the order 
of the universal expansion rate (which comes in through the conversion factor a when 
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we switch the halo profile from physical to comoving coordinates). Thus the condition 
()6.10|) is roughly equivalent to 

TTl I / TTl \ 

k > k^Us) ~ Y^^a'His) - 2y'afi„,o + a^Afi h Mpc~\ (6.15) 

where H{s) is the Hubble expansion parameter at time s. Since fcmin ~ ^fs at most 
times, we see that equation ()6.12p is indeed applicable to all k modes larger than the 
free-streaming wave vector kfs- 

Unfortunately, the opposite k <^ kfs limit has no simple approximate solution 
because of the complicated dependence of the scale factor a on the new time variable s. 
However, we find the following formula to give a decent fit to the solution of ()6.9j) for a 
wide range of k, 

Uk, s) ^ (fe^ /r + fc)2 ^-(^' ^) = ^(^is'k, s)Lik, s), (6.16) 



with 



m,0 



6m{k, s) 



ds'a{s')6m{k, s'){s — s' 



(6.17) 



Typically, F ~ 1, such that for k <C kfs, the growth of 6i, approximately matches that 
of 6m- Therefore, equation ()6.16|) is roughly equivalent to 

p^{x) ~ klK{kfsX)-k p^{x), (6.18) 

in real space, with K{x) = J^^^[K{k)] acting like a normalised filter function with 
window width k^^, that gently smears out the neutrino density contrasts on scales 
below ~ kg^ relative to their CDM counterparts. We shall be using equation ()6.18p 
again in section 18.21 

6.2. Comparison with N-l-body results: when and how the linear theory fails 

Comparing the linear results from this section with A^-l-body simulations from sectional 
it is immediately clear in Figure ^ that the former systematically underestimates the 
neutrino overdensities over the whole range of neutrino and halo masses considered 
in this study. The discrepancy is most prominent in the dense, inner regions 
(r < 1 h~^Mpc), and worsens as we increase (i) the neutrino mass m^, and (ii) the 
halo mass Mvir. The worst case corresponds to when both ruu and Myir are large; the 
case of {m^ = 0.6 eV, M^ir = 0.7x 10^^ H'^Mq}, for example, sees the A^-l-body and the 
linear overdensities differ by a factor of about six. For smaller neutrino and halo masses, 
concordance between the two approaches improves as we move to larger radii. Indeed, 
for {m,^ = 0.15 eV, Mvir = 0.7 x 10^^ /i^^Mq}, complete agreement is seen throughout 
the region of interest. Upon closer inspection, one finds that the linear theory ceases 
to be a faithful approximation once the neutrino overdensity reaches a value of about 
three or four. This is of course fully consistent with the standard lore that perturbative 
methods fail once the perturbations exceed unity and nonlinear effects set in. 
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Can the linear approximation be salvaged? Yes, provided we impose a great deal 
of smoothing. In the case of {m^, = 0.6 eV, Mvir = 0.7 x 10^^ H'^Mq}, for example, 
the overdensities computed from the two different approaches can be reconciled with 
each other if we smooth them both over a scale of roughly 5 h~^Mpc. Such a large 
smoothing length will render the linear method completely useless for the study of 
neutrino clustering on sub-galactic scales (unless of course the neutrino mass is so small 
that the overdensity does not exceed unity by much anyway). But the method can 
still be useful for obtaining quick estimates of Uy/fiy on larger scales for absorption and 
emission spectroscopy calculations (cf. section ing . 

Finally, we note that the neutrino overdensities in Figure 2 of SM are at odds 
with our linear results in Figure H This discrepancy cannot be ameliorated by simply 
supposing that SM have normalised their neutrino densities for three flavours to the 
one flavour average n,^ = np ~ 56 cm~'^, since this normalisation will render some of 
their results — specifically, where riy/ny < 3 — unphysical. Only a normalisation to three 
flavours can give these results a physical meaning, but at the expense of incompatibility 
with our Figure [TJ as well as with SM's own Figure 3. We therefore conclude that SM's 
results as presented in their Figure 2 are erroneous. 

7. Relic neutrinos in the Milky Way 

In this section, we consider relic neutrino clustering in the Milky Way. We compute 
explicitly the number of neutrinos and their distribution in momentum space in our 
local neighbourhood at Earth (r^ ~ 8 kpc from the Galactic Centre). This information 
is essential for any direct search experiment. 

7.1. Background, basic set-up and assumptions 

We perform this calculation using the A^-l-body method of section El but with a few 
modiflcations to the external potential. Firstly, we note that the central region of the 
Milky Way ( < 10 kpc) is dominated by baryonic matter in the form of a disk, a bulge, 
and possibly a rapidly rotating bar [77j. Each of these components has its own distinct 
density profile. Furthermore, in the standard theory of hierarchical galaxy formation, 
baryons and dark matter are initially well mixed, and collapse together to form halos via 
gravitational instability. Galactic structures arise when the baryons cool and fall out of 
the original halo towards the centre [ZH]- As the baryons condense, their gravitational 
forces tend to pull the dark matter inward, thereby distorting the inner ~ 10 kpc of 
the original halo profile (e.g., [ZHl IHUl IMl IH2|)- This kind of modification to the mass 
distribution is important for us at r^^. Fortunately for our calculation, gravitating 
neutrinos do not distinguish between halo and baryons. Therefore, it suffices to use 
simply the total mass distribution inferred from observational data (e.g., rotation curves, 
satellite kinematics, etc.), without any detailed modelling of the individual components. 
What is still missing, however, is the redshift dependence of the mass distribution. 
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Unfortunately, we have not been able to find in the literature any simple parametric form 
for this dependence. However, mass modelling of the Milky Way jH^l Ell suggests that 
certain observationally acceptable bulge+disk+halo models are indeed consistent with 
the aforementioned theory of baryonic compression, and can be traced back to halos 
originally of the NFW form by supposing that the system has undergone a phase of 
adiabatic contraction [7^]. Thus, for our investigation, it is probably fair to think of the 
NFW profile as the initial mass distribution, and the evolution as a smooth transition 
from this initial distribution to the present day one. 

Instead of modelling this transition, however, our strategy here is to conduct two 
series of simulations, one for the present day mass distribution of the Milky Way 
(MWnow) which we assume to be static, and one for the NFW halo (NFWhalo) that 
would have been there, had baryon compression not taken place. The real neutrino 
overdensity should then lie somewhere between these two extremes. 

For the NFWhalo run, we use the parameters Myir = 1 x lO^^M© and c = 12. These 
numbers are taken from the paper of Klypin, Zhao and Somerville (hereafter, KZS) |84j . 
from their "favoured" mass model of the Milky Way. Note that we are not using the 
c-Mvir relation ()4.5|1 . which, as we discussed before, is only a statistical trend. However, 
the concentration parameter c should still carry a redshift dependence a la equation 
()4.6|) . in order to reproduce the correct time dependence of the density profile. 

For MWnow, we adopt the total mass distribution (halo+disk+bulge) presented in 
Figure 3 of KZS, and fit it approximately to a power law from r = to 20 kpc. 



where M(r) means the total mass contained within a radius r. We assume this fit to hold 
for the region inside a physical radius of 20 kpc at all times, i.e., Mfit(r, z) = Mfit(ar, z = 
0). The region outside of this 20 kpc sphere is not affected by baryonic compression 
according to the KZS mass model (cf. their Figure 7) so we adopt the original NFW 
density profile outwardly from 20 kpc. Thus, schematically, we have 



where tq = 20 kpc, Mnfw(^) is the mass contained in an NFW halo at radius r [or 
M\isio{r) in equation ()A.7|l ]. and Mfit(ro) ~ 2 x MnfwI'^o) for the parameters used in 
this analysis. 

7.2. Results and discussions 

Our Milky Way simulation results for four neutrino mass m^, = 0.15, 0.3, 0.45, 0.6 eV are 
displayed in Figure El The shaded region in each plot corresponds to a possible range 
of overdensities at z = 0. At first glance, it may seem unphysical that the apparently 
static MWnow potential (in physical coordinates) should capture so many neutrinos. To 
resolve this "paradox", one must remember that neutrino clustering is studied in the 




1.19367 



M. 



(7.1) 



M(r,r <ro) = Mfit(r), 

M(r, r > ro) = M^F^{r) - MNFw(ro) + Mfit(ro), 



(7.2) 
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Figure 5. Relic neutrino number density per flavour, n,y = np, in the Milky Way for 
various neutrino masses. All curves are normalised to n,^ = fip ~ 56 cm""^. The top 
curve in each plot corresponds to the MWnow run, and the bottom to the NFWhalo 
run. The enclosed region represents a possible range of overdensities at z = 0. 



context of an expanding universe; the (unbound) neutrino thermal velocity decreases 
with time [equation ()4.H) ]. thus causing them to be more readily captured. Equivalently, 
in comoving coordinates, it is easy to see that while the neutrino conjugate momentum 
()3.2|) does not redshift, the MWnow potential well shrinks in size and deepens with time. 

In each scenario we studied, the final momentum distribution at is almost 
isotropic, with a zero mean radial velocity (vr), and second velocity moments that 
satisfy approximately the relation 2(f^) = {v^) (cf. Table Ej). For this reason, we plot 
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Table 2. Velocity moments at for various neutrino masses in the MWnow and 
NFWhalo runs (see text for definitions). The first column shows the overdensities 
nu/fiu- The second, third and fourth columns show the mean radial, transverse and 
absolute velocities in terms of the diniensionless quantities (j/r), {vt) and (y), where 
y — m^v /Tjjfl. In the last three columns are the second moments. The corresponding 
values for a relativistic Fermi-Dirac distribution are displayed in the first row. 







(2/r) 


(2/t) 


(y> 








Relativistic Fermi-Dirac 


1 





2.48 


3.15 


4.31 


8.63 


12.94 


NFWhalo, 


= 0.6 eV 


12 


0.0 


3.4 


4.3 


6.9 


13 


20 


NFWhalo, 


= 0.45 eV 


6.4 


0.0 


2.8 


3.5 


4.6 


9.5 


14 


NFWhalo, 


= 0.3 eV 


3.1 


0.0 


2.3 


3.0 


3.6 


7.3 


11 


NFWhalo, 


= 0.15 eV 


1.4 


0.0 


2.3 


2.0 


3.8 


7.6 


11 


MWnow, = 


= 0.6 eV 


20 


0.0 


4.0 


5.1 


9.3 


18 


28 


MWnow, niy = 


= 0.45 eV 


10 


0.0 


3.1 


4.0 


6.1 


12 


18 


MWnow, 771^ = 


= 0.3 eV 


4.4 


0.0 


2.5 


3.2 


3.9 


8.0 


12 


MWnow, nil, = 


= 0.15 eV 


1.6 


0.0 


2.3 


2.9 


3.7 


7.3 


11 



the smoothed, or coarse-grained, phase space densities f{r^,p) only as functions of the 
absolute velocity (cf. Figure IHI). 

As expected, the coarse-grained distribution fir^.p) for the case with the highest 
overdensity (MWnow, rui, = 0.6 eV) resembles the original Fermi-Dirac spectrum the 
least, while / for the case with the lowest overdensity (NFWhalo, m^, = 0.15 eV) 
is almost Fermi-Dirac-like. All spectra share the feature that they are flat at low 
momenta, with a common value of ~ 1/2. The turning point for each distribution 
coincides approximately with the "escape momentum" pesc (i-e., times the escape 
velocity fesc = ^2|0(r0)|) for the system concerned, beyond which the phase space 
density falls off rapidly, until it matches again the Fermi-Dirac function at the very high 
momentum end of the spectrum. Deviation from the original Fermi-Dirac spectrum is 
therefore most severe around pesc- 

The maximum value of / is a little less than 1/2. This is consistent with the 
requirement that the final coarse-grained density must not exceed the maximal value 
of the initial fine-grained distribution, / < max(/o) jHS]- For neutrinos, /o has a value 
of 1/2 at p = 0. Thus, our / not only satisfies but completely saturates the bound at 
low momenta up to Pcsc? forming a kind of semi-degenerate state that can only be made 
denser by filling in states above Pesc-% However, since neutrinos with momenta above pesc 
do not become gravitationally bound to the galaxy/halo, these high momentum states 
are much less likely to be fully occupied. This explains /'s rapid drop beyond Pesc- Also, 
the hottest neutrinos are not significantly affected by the galaxy/halo's gravitational 
forces. Therefore the very high end of the momentum distribution remains more or less 
Fermi-Dirac-like. Finally, we note that because the filling of phase space happens from 

% This degeneracy should not be confused with that arising from the Pauli exclusion principle. 
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Figure 6. Momentum distribution of relic neutrinos at for various neutrino masses. 
The red (solid) line denotes the MWnow run, while the dashed line represents the 
NFWhalo rim. The relativistic Fermi-Dirac function is indicated by the dotted line. 
The escape velocity Wosc = \/2[^(r®)| is 490 km s~^ and 450 km s^^ for MWnow 
and NFWhalo respectively, corresponding to "escape momenta" j/csc = fn^Vcsc/T^,o of 
(5.9, 4.4, 3.0, 1.5) and (5.4, 4.1, 2.7, 1.4) for = (0.6, 0.45, 0.3, 0.15) eV. 
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bottom up, the mean momenta for the least clustered cases tend to be lower than the 
Fermi-Dirac value (p) ^ 3.15T,^_0) in contrast with the naive expectation that clustering 
is necessarily accompanied by an increase in the neutrinos' average kinetic energy. 



7.3. Tremaine-Gunn hound 

It is interesting to compare our results with nominal bounds from phase space arguments. 
By demanding the final coarse-grained distribution to be always less than the maximum 
of the original fine-grained distribution, Tremaine and Gunn |86j argued in 1979 that if 
neutrinos alone are to constitute the dark matter of a galactic halo, their mass must be 
larger than 20 eV, assuming that the halo has a Maxwellian phase space distribution 
as motivated by the theory of violent relaxation [HHl IHH IHH]- A modern version of 
this bound, in which the assumption about the phase space distribution is relaxed and 
which allows for contribution to the total gravitational potential from more than one 
form of matter, can be found in reference IHS] ■ The revised mass bound may be written, 
alternatively, in the form of a constraint on the overdensity, which reads 

3 3 

n. 9C(3)T3o- 

For the neutrino masses m^, = (0.6,0.45,0.3,0.15) eV, this expression evaluates to 
(19,8.0,2.4,0.3) and (15,6.4,1.8,0.25) for MWnow and NFWhalo, respectively, at r®. 

At first sight, some of our numerical results seem to have completely violated the 
bound ()7.3p . But this cannot be, since we have seen explicitly that all of our final 
coarse-grained distributions satisfy perfectly the constraint / < max(/o). Furthermore, 
an upper bound of 0.25 on the overdensity for a 0.15 eV neutrino is obviously nonsensical. 

The answer, as the astute reader would have figured, lies in accounting. In the 
derivation of ()7.3|1 . a semi-degenerate distribution has been summed only up to the 
momentum state corresponding to the escape velocity of the system. Neutrinos with 
higher momenta that could be hovering around in the vicinity have been completely 
ignored. In contrast, in our calculations, it is of no concern to us whether or not the 
relic neutrinos actually form bound states with the galaxy/halo. Therefore it is more 
appropriate for us to sum every neutrino in sight, rather than imposing a cut-off at 
Vesc- However, if we had imposed such a cut-off, one can easily see from Figure IHl that 
our overdensities would have just saturated the bound ()7.3j) . so there is no conflict. 
Nonetheless, this illustrates how nominal bounds such as ()7.3p must be used with care. 

Before we conclude this section, let us note that for an NFW halo, the bound ()7.3j] 
can be written as 



3 



< 



^ 9C(3)T3 



2GM^i, ln(l + r/r. 



-1 3/2 



- 0.58 X — ^ 




M^ir \ fh-'Mpc\ ln(l + r/r,)]'/' 



,eVy [{lO^'^h'mQj \ r J g{c) 
where the function g{c) is defined in equation ()A.8|) in the Appendix. Figure [7| shows 
the bound as a function of radius for the various halo and neutrino masses considered 
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Figure 7. The Tremaine-Gunn bound on the neutrino overdensity for various halo 
and neutrino masses (dashed hnes). The red (sohd) hnes correspond to our A^-l-body 
results from section IHl 

in section along with the overdensities obtained from A^-l-body simulations. We 
find the limit ()7.3p to be saturated for the lightest halo and neutrino masses (and also 
some apparent violation of the bound due to different accounting). This explains why 
some of the overdensities exhibit an almost ml dependence (cf. Figure |3I). On the 
other hand, the heaviest {rui,, Mvir} set is short of the bound by at least one order of 
magnitude. This is also consistent with expectations: heavier halo and neutrino masses 
give a higher escape momentum pesc, and higher momentum states are more difficult 
to fill up to a semi-degenerate level, since there are less particles in these states in the 
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original Fermi-Dirac distribution to begin with. 



8. Implications for detection 

In this section, we determine the implications of our clustering results for the direct 
detection of relic neutrinos, in contrast with the purely cosmological inferences discussed 
in section El We consider various proposed detection methods based on scattering 
processes, involving the relic neutrinos both as a beam and as a target. In particular, 
we shall discuss (i) coherent elastic scattering of the relic neutrino flux off target 
matter in a terrestrial detector f section l8.1|) . as well as (ii) the scattering of extremely 
energetic particles (accelerator beams or cosmic rays) off the relic neutrinos as a target 
(sectioning). 



8.1. Flux detection 

The low average momentum (p) = {y) T^^q of the relic neutrinos (cf. Table 12]) 
corresponds to a (reduced) de Broglie wavelength of macroscopic dimension, X = l/{p) = 
0.12 cm./ {y) (cf Table Ej). Therefore, one may envisage scattering processes in which 
many target atoms act coherently |31 |3] over a macroscopic volume A^^, so that the 
reaction rate for elastic scattering becomes proportional to the square of the number of 
target atoms in that volume. Compared to the case where the neutrinos are elastically 
scattered coherently only on the individual nuclei of the target, the rate in this case is 
enhanced by a huge factor of 

where Na is the Avogadro constant, A is the atomic mass, and pt is the mass density 
of the target material."'" 

By exploiting the above coherence effect, a practical detection scheme for the local 
relic neutrino flux is based on the fact that a test body of density pt at Earth will 
experience a neutrino wind force through random neutrino scattering events, leading to 
an acceleration [U El IHl 1^ 

47r 

'^'^ flux mom. transfer 

.2X10- (!^) (^) (^) (^77^)' (8.2) 

\n^J \ Vrc\ ) \g/cm^y \h/[m^v^c\) J 

where a^N — G"]^ fn^/T^ is the elastic neutrino-nucleon cross section, and frci = (1*^ — i^el) 
is the mean velocity of the relic neutrinos in the rest system of the detector. Here, 
v@ ~ 2.3 X 10^ km s~^ ~ 7.7 x 10~'^c denotes the velocity of the Earth through the 

In the case of coherent scattering, it is possible, in principle, to measure also the scattering amplitude 
itself |9U[l91ll^ . which is linear in the Fermi coupling constant Gp- However, one needs a large lepton 
asymmetry for a non-negligible effect. 
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Table 3. Properties of the rehc neutrinos at for various neutrino masses from 
the MWnow and NFWhalo runs (see text for definitions) relevant for their direct 
detection. The first column shows the overdensities. Columns two, three and four show, 
respectively, the mean absolute momenta, the associated mean reduced de Broglie 
wavelengths and the mean absolute velocities (in units of c) . The corresponding values 
for a relativistic Fermi-Dirac distribution are displayed in row one. 







(P> 




{v) 


Rclativistic Fermi-Dirac 


1 


5.3 X 10-4 eV 


3.7 X 10-2 cm 


see (|4.1f) 


NFWhalo, TO^ = 0.6 eV 
NFWhalo, = 0.45 eV 
NFWhalo, = 0.3 eV 
NFWhalo, = 0.15 eV 


12 
6.4 
3.1 
1.4 


7.2 X 10-4 eV 
5.9 X 10-4 eV 
5.0 X 10-4 eV 
3.4 X 10-4 eV 


2.7 X 10-2 cm 
3.4 X 10-2 cm 
3.9 X 10-2 cm 
5.9 X 10-2 cm 


1.2 X 10-3 

1.3 X 10-3 
1.7 X 10-3 
2.2 X 10-3 


MWnow, m„ = 0.6 eV 
MWnow, = 0.45 eV 
MWnow, = 0.3 eV 
MWnow, = 0.15 eV 


20 
10 
4.4 
1.6 


8.5 X 10-4 eV 
6.7 X 10-4 eV 
5.4 X 10-4 eV 
4.9 X 10-4 eV 


2.3 X 10-2 cm 
2.9 X 10-2 cm 
3.7 X 10-2 cm 
4.1 X 10-2 cm 


1.4 X 10-3 

1.5 X 10-3 
1.8 X 10-3 
3.2 X 10-3 



Milky Way. Expression ()8.2j) is valid as long as the radius of the target is smaller than 
the reduced de Broglie wavelength X = h/ {m^Vy-d) of the relic neutrinos. Furthermore, it 
applies only to Dirac neutrinos. For Majorana neutrinos, the acceleration is suppressed, 
in comparison with ()8.2|1 . by a factor of (frei/c)^ ~ 10"^ for an unpolarised target, or 
Vj-ei/c ~ 10"^ for a polarised one. A target size much larger than X can be exploited, 
while avoiding destructive interference, by using foam-like ^ or laminated [5] materials. 
Alternatively, grains of size ~ X could be randomly embedded (with spacing ~ A^) in a 
low density host material jHH IHS] • 

To digest these estimates, we note that the smallest measurable acceleration at 
present is > 10"^^ cm s"^, using conventional Cavendish-type torsion balances. Possible 
improvements with currently available technology to a sensitivity of > IQ-^'^ cm s"^ 
have been proposed jUHlinZI- However, such a sensitivity is still off the prediction ()8.2p 
by at least three orders of magnitude, as an inspection of the currently allowed range 
of local relic neutrino overdensities in Table El reveals. Therefore, we conclude that 
an observation of this effect will not be possible within the upcoming decade, but can 
still be envisaged in the foreseeable future (thirty to forty years according to reference 
P3] . exploiting advances in nanotechnology) , as long as our known light neutrinos are 
Dirac particles. Should they turn out, in the meantime, to be Majorana particles, flux 
detection via mechanical forces will be a real challenge. 

Let us note finally that the background contribution to the acceleration ()8.2|] 
from the solar pp neutrinos [flux ~ 10^^ cm-^s"^, {E,y) ~ 0.3 MeV (e.g., [HH])], 
^i/sun ^ lO"^"^ cm s"^ [0], may be rejected by directionality. The background from 
weakly interacting massive particles (WIMPs with mass m^) jB], 

^wiMP ^ n^^y^ N^A a^N 2 v^ei (8.3) 

flux mom. transfer 
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Table 4. Beam parameters of forthcoming accelerators and expected interaction rates 
with rehc neutrinos. 



accelerator 


N 


A 


Z 


En [TeV] 


L [km] 


/ [A] 


RvA 


LHC 


P 

Pb 


1 

208 


1 
82 


7 

574 


26.7 
26.7 


5.8 X 10-1 
6.1 X 10-3 


2 X 10-8 yr-i 
1 X 10-5 yr-i 


VLHC 


P 

Pb 


1 

208 


1 
82 


87.5 
7280 


233 
233 


5.7 X 10-2 
5.7 X 10-'' 


2 X lO-'^ yr-i 
1 X IQ-'^ yr-i 


ULHC 


P 


1 


1 


10^ 


4 X 10^ 


1.0 X 10-1 


10 yr-i 



-^""^^ ''(o.3GeV/cm3)(lO-3c) ( lOo) ( IQ-^w) "^"^ ' 

should they be the main constituent of the galactic dark matter with mass density 
= n^m^ ~ 0.3 GeV cm~^ at r^, can be neglected as soon as the WIMP-nucleon 
cross section a^N is smaller than ~ 3 x lO—^^ cm^. This should be well established 
by the time relic neutrino direct detection becomes a reality. Note also that neutrinos 
produced in the Earth's atmosphere do not contribute appreciably to at because of their 
small flux; for ~ 0.1 ^ 10 GeV, the integrated flux is < 1 cm'^s"^ (e.g., P^). 

8.2. Target detection 

Let us consider next the idea to take advantage of the fact that, for center-of-mass 
(cm.) energies below the W- and Z-resonances, the weak interaction cross sections 
grow rapidly with energy. One may then contemplate the possibility to exploit a flux of 
extremely energetic particles — either from accelerator beams or from cosmic rays — for 
scattering on the relic neutrinos. 

8.2.1. At accelerators We start with a discussion of the prospects to detect interactions 
with the relic neutrinos at forthcoming accelerator beams such as the LHC |100j and 
the VLHC jlOlj . with beam energies E^i ranging from 7 TeV for the LHC running with 
protons, up to several thousands of TeV for the heavy ion option at the VLHC (cf. 
Table llj. At these beams, the attainable momentum transfers and cm. energies, 

v/i = .J^I^, ^ 4,5 (^) (^) MeV. (8,4) 

are so small that the cross sections for their interactions with the relic neutrinos are 
enhanced by a factor ~ due to coherent elastic scattering over the size of the 
nucleus pjJ2i j. and grow linearly with the beam energy, 

^= Gl./'. - 3.4 X 10-« (;;^) (^) cm^ (8,5) 
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This leads to a scattering rate \WS\ llU4t ll(J5j 

K^N - E ^^^v^TV^ I/iZe) (8.6) 

for a beam of particles ^A^, with charge Z e, length L and current /. In view of the 
currently allowed range of local relic neutrino overdensities displayed in Table 01 and 
the beam parameters of the next generation of accelerators summarised in Table IH 
the expected rate ()8.6p is clearly too small to give rise to an observable effect in 
the foreseeable future. The extremely energetic 574 TeV lead beam at the LHC will 
have less than 10~^ interactions per year with relic neutrinos, for m^, < 0.6 eV and, 
correspondingly, Uv/fiy < 20. Even a lead acceleration option for the VLHC, with 
En — 7280 TeV (or, equivalently, 35 TeV per nucleon) and a current / ~ 5.7 x 10^^ A (a 
hundredth of the nominal current of the p running mode) will give less than 10~^ events 
per year for the most optimistic neutrino mass scenario. Thus, there is little hope, in 
the foreseeable future, to detect relic neutrino using terrestrial accelerator beams. 

Let us nevertheless dream about the far future, in which an Ultimate Large Hadron 
Collider (ULHC) exists and is able to accelerate protons to energies above 10'' TeV* 
in a ring of ultimate circumference L ~ 4 x 10^ km around the Earth, thus leading 
to an interaction rate of more than one event per year (cf. Table E}. Even under 
these most optimistic circumstances, is it possible to reliably detect these interactions? 
Clearly, elastic scattering of the beam particles with the relic neutrinos — one of the 
contributions to the rate ()8.6p — will be next to impossible to detect because of the 
small momentum transfers involved (~ 1 GeV at En ~ 10^ TeV). A very promising 
alternative is to consider again a heavy ion beam, and to exploit the contribution of the 
inverse beta decay reaction, 

iN + u,-.i^,N + e- , (8.7) 

to the rate ()8.6|) . This reaction changes the charge of the nucleus, causing it to follow an 
extraordinary trajectory and finally to exit the machine such that it becomes susceptible 
to detection |104^ll()9] . A detection of this reaction would also clearly demonstrate that 
a neutrino was involved in the scattering. 

8.2.2. With cosmic rays In the meantime, until the ULHC has been constructed, target 
detection of the relic neutrinos has to rely on extremely energetic cosmic rays. In 
fact, cosmic rays with up energies up to E^ ~ 10^° eV have been seen by air shower 

* Note that a collider at this energy has to be built anyhow if one wishes to explore the "intermediate" 
scale (Mew -Mqut)^^^ ^ 10^" GeV between the electroweak scale Mew ^ 1 TeV and the scale of grand 
unification Mqut ~ 10^^ GeV. The intermediate scale is exploited in many schemes of supersymmetry 
breaking (e.g., jlOfij ) and in seesaw mechanisms for neutrino masses |1()71 11(18] . 
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observatories. The corresponding cm. energies are 

yi= V^^=. 14 (^^Y' {j^f GeV (8.8) 

when scattering off the relic neutrinos. These energies are not too far from the W- and 
Z-resonances, at which the electroweak cross sections become sizeable. Indeed, it was 
pointed out long ago by Weiler [71IH] (for earher suggestions, see II 111 11121111^ II 14j ) 
that the resonant annihilation of extremely energetic cosmic neutrinos (EECi/) with relic 
anti-neutrinos (and vice versa) into Z-bosons appears to be a unique process having 
sensitivity to the relic neutrinos. On resonance, 

^res ^ il4 ^ 4 ^ ^q21 f ^\ ^ (g 

the associated cross section is enhanced by several orders of magnitude, 

((Xann) = J ds / m\ (^uM ^ 2'kV2Gf ^ 4 X 10"^^ cm^ (8.10) 

leading to a "short" mean free path 1^ = {n,, (cTann))^^ — 1.4 x 10^ Mpc which is 
only about 48 h times the Hubble distance. Neglecting cosmic evolution effects, this 
corresponds to an annihilation probability for EECz/ from cosmological distances on the 
relic neutrinos of 2h~^%. 

The signatures of annihilation are (i) absorption dips [3 El El (see also |115| lll(i| 
I117j ) in the EECz/ spectrum at the resonant energies, and (ii) emission features |10| 
ITT| IT2l IT3l IT^ (Z-bursts) as protons (or photons) with energies spanning a decade 
or more above the predicted Greisen-Zatsepin-Kuzmin (GZK) cutoff at Eqzk — 
4 X 10^^ eV |118t I119j . This is the energy beyond which the CMB is absorbing to 
nucleons due to resonant photopion production. jj 

The possibility to confirm the existence of relic neutrinos within the next 
decade from a measurement of the aforementioned absorption dips in the EECi/ 
flux was recently investigated in j9j. Presently planned neutrino detectors (Pierre 
Auger Observatory IceCube [HS], ANITA [121], EUSO [123, OWL |I2H1, and 

SalSA |129j ) operating in the energy regime above 10^^ eV appear to be sensitive 
enough to lead us, within the next decade, into an era of relic neutrino absorption 
spectroscopy, provided that the flux of the EECz/ at the resonant energies is close to 
current observational bounds and the neutrino mass is sufficiently large, > 0.1 eV. 
In this case, the associated Z-bursts must also be seen as post-GZK events at the 
planned cosmic ray detectors (Pierre Auger Observatory, EUSO, and OWL). 

What are the implications of relic neutrino clustering for absorption and emission 
spectroscopy? Firstly, absorption spectroscopy is predominantly sensitive to the relic 
neutrino background at early times, with the depths of the absorption dips determined 
largely by the higher number densities at large redshifts {z ^ 1). Since neutrinos do 
not cluster significantly until after z < 2, clustering at recent times can only show up as 

The association of Z-bursts with the mysterious cosmic rays observed above Eqzk is a controversial 

possibility [rmiminiorra nnniiTTnimiT^ . 
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Figure 8. "Large scale" overdensities (« = ly, CDM) in tlie local universe, with the 
Milky Way at the origin, estimated from equation (|6.18|l assuming the function K{x) 
to be a Gaussian. The black (solid) line corresponds to the local CDM distribution 
inferred from peculiar velocity measurements |13H (see also |132| ') smeared over the 
surface of a sphere with radius r [21 • The dotted line is the neutrino overdensity for 
m^, = 0.6 eV, short dash 0.3 eV, long dash 0.15 eV, and dot-dash 0.04 eV. 

secondary dips with such minimal widths in energy jl30j that they do not seem hkely 
to be resolved by planned observatories. 

On the other hand, emission spectroscopy is directly sensitive to the relic neutrino 
content of the local universe {z < 0.01 -v^ ?"gzk < 50 Mpc). However, since the neutrino 
density contrasts approximately track those of the underlying CDM above the neutrino 
free-streaming scale k^^ (cf. section IH)), it is clear that there cannot be a substantial 
neutrino overdensity over the whole GZK volume (~ Tgzk)- Indeed, if we take the linear 
fitting formula ()6.18|) . and apply it to the local CDM distribution inferred from peculiar 
velocity measurements (with a ~ 5 Mpc smoothing), one can see in Figure |H1 that the 
estimated neutrino overdensity is always < 2. Hence the overall emission rate cannot 
be significantly enhanced by gravitational clustering. 

Or we could imagine doing "relic neutrino tomography" of the local universe, by 
exploiting the fact that there are several galaxy clusters (> lO^^M©), such as Virgo 
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(distance ~ 15 Mpc) and Centaurus (~ 45 Mpc), within the GZK zone with significant 
neutrino clustering (cf. Figure Q). One could conceivably search for directional 
dependences in the emission events as a signature of EECz/ annihilating on relic anti- 
neutrinos (and vice versa). For example, AGASA has an angular resolution of ~ 2° 
|133j . This is already sufficient to resolve the internal structures of, say, the Virgo cluster 
(distance ~ 15 Mpc, Mvir ~ 8 x IO^^Mq) which spans some 10° across the sky. Using our 
A^-l-body clustering results in Figure ^ the average neutrino overdensity along the line 
of sight towards and up to Virgo is estimated to be ~ 45 and ~ 5 for nii, = 0.6 eV and 
0.15 eV respectively, given an angular resolution of ~ 2°. The corresponding increases 
in the number of events coming from the direction of the Virgo cluster relative to the 
unclustered case, assuming an isotropic distribution of EECz/ sources, are given roughly 
by the same numbers, since protons originating from ~ 15 Mpc away arrive at Earth 
approximately unattenuated. The numbers improve to ~ 55 and ~ 8 respectively with 
a finer ~ 1° angular resolution. 

Note that our estimates here are generally a factor of few higher than the predictions 
of SM. This is expected, because the linear method adopted in their analysis cannot 
account for additional clustering from nonlinear effects, as we demonstrated in section IHl 

9. Conclusion 

We have conducted a systematic and exhaustive study of the gravitational clustering 
of big bang relic neutrinos onto existing CDM and baryonic structures within the fiat 
ACDM model, with the aim of understanding their clustering properties on galactic 
and sub-galactic scales for the purpose of designing possible scattering-based detection 
methods. Our main computational tools are (i) a restricted, A^-l-body method 
(section E)), in which we neglect the gravitational interaction between the neutrinos 
and treat them as test particles moving in an external potential generated by the 
CDM/baryonic structures, and (ii) a semi-analytical, linear technique (section Ej), which 
requires additional assumptions about the neutrino phase space distribution. In both 
cases, the CDM/baryonic gravitational potentials are calculated from parametric halo 
density profiles from high resolution A^-body studies jHOl IS] (section and/or from 
realistic mass distributions reconstructed from observational data (e.g., [5^ I131j ). 

Using these two computational techniques, we track the relic neutrinos' accretion 
onto CDM halos ranging from the galaxy to the galaxy cluster variety {M^ij. ~ 10^^ — >■ 
10^^ Mq), and determine the neutrino number densities on scales ~ 1 ^ 1000 kpc for 
neutrino masses satisfying current constraints 1)2.11) from CMB and LSS (Figures^ and 
12}. Because we can simulate only a finite set of halo and neutrino parameters, we 
provide also additional plots illustrating the approximate dependences of the neutrino 
overdensities on the halo and neutrino masses (Figures El and 0}. These can be used for 
interpolation between simulation results. Furthermore, we find that the linear technique 
systematically underestimates the neutrino overdensities over the whole range of halo 
and neutrino masses considered in this study (Figure Q). Reconciliation with A^-l-body 
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simulations can only be achieved if we impose a smoothing scale of > 1 Mpc, or if the 
overdensity is no more than three or four. We therefore conclude that the linear theory 
does not generally constitute a faithful approximation to the Vlasov equation in the 
study of neutrino clustering on galactic and sub-galactic scales ( < 50 kpc). However, 
it may still be useful for finding the minimum effects of neutrino clustering in other 
contexts not considered in this work (e.g., the nonlinear matter power spectrum |134j ). 

Next we apply our A^-l-body method to calculate the relic neutrino number 
density in the Milky Way (section Ej), especially their phase space distribution in our 
local neighbourhood at Earth r^, taking also into account contributions to the total 
gravitational potential from the galactic bulge and disk. We find a maximum overdensity 
of ~ 20 per neutrino flavour in our immediate vicinity, provided that the neutrino mass is 
at its current upper limit of 0.6 eV (Tableland Figure Ej). For neutrino masses less than 
0.15 eV, the expected overdensity from gravitational clustering is less than two. The 
associated coarse-grained momentum spectra show varying degrees of deviation from 
the relativistic Fermi-Dirac function, but share a common feature that they are semi- 
degenerate, with phase space density / ~ 1/2, up to the momentum state corresponding 
to the escape velocity from the Milky Way at (FigureE))- This means that the neutrino 
number densities we have calculated here for are already the highest possible, given 
the neutrino mass, without violating phase space constraints (e.g., [S3 EH])- In order to 
attain even higher densities, one must now appeal to non-standard theories (e.g., |135j ). 

In terms of scattering-based detection possibilities, this meager enhancement in 
the neutrino number density in the Milky Way from gravitational clustering means 
that relic neutrinos are still far from being detected in fully earthbound laboratory 
experiments. For flux detection methods based on coherent elastic scattering of relic 
neutrinos off target matter in a terrestrial detector f section iS.l\\ . a positive detection 
could be thirty to forty years away [95^ , provided that light neutrinos are Dirac particles. 
For light Majorana neutrinos, another ~ 10'^ times more sensitivity would be required 
in the detector for a positive signal. Target detection methods using accelerator beams 
(section 18.2.11) seem equally hopeless, unless the accelerator is the size of the Earth and 
operates at an energy of ~ 10^ TeV (Table EJ. 

Meanwhile, target detection using extremely energetic cosmic neutrinos (EECz/, 
> 10^^ eV) remains the only viable means to confirm the existence of big bang relic 
neutrinos within the next decade or so f section l8.2.2j) . Resonant annihilation of EECz/ 
on relic neutrinos can be revealed as absorption dips in the EECz/ flux (e.g., PJ, or as 
emission features in the Z-decay products. However, since absorption spectroscopy is 
largely insensitive to late time {z < 2) relic neutrino clustering, our findings here have 
little impact on the conclusions of 9 . On the other hand, emission spectroscopy is 
sensitive to the relic neutrino content of the local GZK zone, Vgzk ~ 50^ Mpc^. While 
we find no significant large scale clustering within Vgzk (Figure |H1) and therefore no 
significant enhancement in the overall emission rates, it is still conceivable to exploit the 
considerable neutrino overdensities in nearby galaxy clusters, and search for directional 
dependences in the post-GZK emission events. For the Virgo cluster, for example, 
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we estimate the event rate from the central 1° region to be ~ 55 and ~ 8 times the 
unclustered rate for neutrino mass m,^ = 0.6 eV and 0.15 eV respectively, assuming 
an isotropic distribution of EECz/ sources. [Our estimates differ by a factor of a few 
from the predictions of because their linear technique fails to account for additional 
clustering from nonlinear effects, as we have demonstrated in this study (section IH))]. 
Planned observatories such as the Pierre Auger Observatory |124j . EUSO |127j and 
OWL |128j will have sufficient angular resolution to, in principle, see this enhancement. 
However, considering the rapidly improving constraints on both the EECz/ flux and 
neutrino masses, it remains to be seen if the enhancement can indeed be observed with 
enough statistical significance. 
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Appendix A. Details of A^-l-body simulations 

We give in this appendix a detailed discussion of the techniques used in our A^-l-body 
simulations. 



Appendix A. 1. Hamilton's equations 



Following inn], we begin with a Lagrangian L{x, x,t) = a \myv'^ — mj,0(|a;|, r) for a 
test particle (neutrino) with mass my and peculiar velocity v = x moving in a spherically 
symmetric gravitational potential well 0(a3, r) generated by the CDM halo. Because of 
spherical symmetry, the motion of the test particle is confined to a plane. This allows 
us to switch to polar coordinates {r, 6'} to obtain 

■1 



L(r,^,r,^,r) 



-mJr"^ + r'^9'^) 



m,. 



r, r 



(A.i; 



and eliminate two superfiuous variables in the process. The canonical momenta 
conjugate to r and 6 are pr = dL/dr = am^r and i = rpg = dL/dO = am^r'^O 
respectively, leading to the Hamiltonian 



H{r,9,pr,i,T) 



1 



2am„ 



and hence Hamilton's equations 

dr Pr 
dr amy ' 

dpr f' _ d(f) 
dr 



de_ 

dr 



+ amy(j){r, r). 



(A.2) 



am,. 



amyT" 



dr 



am,yr 



2' 



(A.3) 



where the last equality expresses the conservation of conjugate comoving angular 
momentum £. 
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For computational purposes, it is convenient to rewrite the set of differential 
equations ()A.3|) in terms of redshift z and the normalised quantities Ur = Pr/^u = 
and U0 = a. /my = ar'^6, thus rendering ()A.3j) into the form 

dr Ur dO Ug 



dz d ^ dz dr^ 

dur ^ f^l 2^^\ ^'^t 
dz d\r^ dr / dz 



dUr lfuj_^,d^\ du^^^^ 



where a = {1 + z)~^ and d = Hq^g'^ fi^.o + flA,o are, respectively, the scale factor 
and the expansion rate, the latter of which is determined by the present day Hubble 
constant Hq and the cosmological model in hand. 

The gravitational potential 0(r, r) obeys the Poisson equation 

- 1|: (r^^) = AnGa%,.ir)6Ur, r). (A.5) 

Since we are treating the CDM halo as a density perturbation sitting on top of a uniform 
background, i.e., Pm{T)Sm{r,T) — > phaio(^), we have the following (partial) solution: 

90 _ AirCa ^ r G 

dr 

where Mhaio('") is the physical mass contained in a sphere of radius r. For an NFW halo, 

Mhaio(r) = Ana^p^r^^gix) = M^i,g{x)/ g{c), x = r/r„ (A.7) 

^(x) = ln(l + x)--^, (A.8) 
1 + X 



r Pi.Mr', r)r"dr' = ^M^^^r), (A.6) 
Jo ar^ 



9 / M^;r \ 



0, 



1 + 2 \1.5X 10l3/i-lM(r 



(A.9) 



1/3 



= 9.51 X 10-1 c-'Q;^fA;.^/' ( iqi^I-imJ ^"'^P"' ^^"^^^ 
The meanings of these various symbols are explained in section |3] 

Appendix A. 2. Discretising the initial phase space distribution 

The number of neutrinos at some initial time in the interval {x,p) {x + dx,p + dp) 
is defined to be 

dN = f{x,p,Ti)d^x d^p. (A.ll) 

We assume the initial phase space distribution to be isotropic and homogeneous, i.e., 
f{x,p,Ti) = foip). This initial homogeneity implies that the ensemble's subsequent 
evolution under a spherically symmetric potential will depend only on the initial radial 



distance r, radial momentum pr, and transverse momentum pt = yp"^ —p"^-- Thus we 
may rewrite the phase space volume element as 

d^x sin 6' dO dcj) dr, d^p — ^ px dpx dp^ dip, (A. 12) 
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with r G [0, oc), 9 e [0,7i], (/) e [0, 2n), and pr G [0, oc), p,. G (-oo, oo), ip G [0, 27r). 

In an ideal and perfectly deterministic calculation, one would track the motion 
of every particle for every occupied combination of r, pr and pt- This is obviously 
impossible in practice; one must therefore resort to simulating only a representative set 
of particles and endow each with a "weight" according to the phase space region from 
which the particle is drawn. We show below how this weight is calculated. ff 

Suppose that a point particle initially at (r, p^^Pt) is representative for all particles 
in the initial phase space interval {ra,Pr,a,PT,a) — ^ {j'b,Pr,b,PT,b)- The weight carried by 
this representative particle is defined as 

Wi= / dN, (A. 13) 

where Jg^^^^ means summing over all 9, and ip, which is simple: Jg^^^^, sin 9 d9 dcj) dip = 
Stt^. The spatial integral over r is also readily calculable. The remaining phase space 
integral Jp^'^^'p^'^ f {p)pTdpTdpr can be solved by way of the parameterisation Pr = p cosip, 
and Pt = psinip, and hence 

Pt dpTdpr p^ simp dip, p G [0, oo), ip G [0, tt). (A. 14) 

Thus the weight carried by a point particle centred on and representing the interval 

ira,Pa,^a) ^ irb,Pb,^b) IS glveu by 

fl'b fVb fPb 

u;, = 87rX'o/ r'dr f{y)y^dy simpdiP, (A.15) 

Jra Jya J^a 

where y = p/Ty^ is a dimensionless momentum variable. 

Appendix A. 3. Kernel method for density profile estimation 

We construct the neutrino density profiles from the discrete outputs of our A^-l-body 
simulations using the kernel method of and jZOj, which we reproduce here for 
completeness. 

The number density corresponding to a set of particles is estimated to be 



1=1 



1 , 
\r 



-.J, (A.16) 

where Wi is the weight carried by the ith particle, h is the window width, and is a 
normalised kernel such as the Gaussian kernel. 

In order to obtain a spherically symmetric profile, we smear every particle around the 
surface of a sphere with radius centred on r = 0. The spherically symmetric density 
estimate is then 

n{r)=J2^K{r,r„h), (A.18) 
i=i 

ff We note for completeness that the concept of super particles with equal masses or, in our language, 
weights is adopted in most full scale A^-body simulations (T5] . The representative point particles are 
drawn randomly from the initial phase space distribution. 
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with 

K{r,ri,h) = ^ j d<p j dO sin OK (^h~^^Jr^ + rf - 2rri cos 6^ . (A. 19) 
Substituting for the Gaussian kernel, we get 

'<ir. n. h) = ^pl,^^ [e-i'-r.rU"' - e-<'-.)V-'] . (A.20) 

The optimal choice of h generally depends on the underlying function n(r) |70j. In our 
analysis, however, a constant h is adopted for simplicity. Our momentum distributions 
in section [7| are also constructed with this kernel method. 
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